v20200408.1936

Objetivos

  • reconhecer e mencionar propriedades da distribuição \(F\).
  • reconhecer as indicações e aplicar ANOVA para três ou mais condições independentes.
  • reconhecer as indicações e aplicar ANOVA para três ou mais condições dependentes (medidas repetidas).
  • definir hipóteses estatísticas nula e alternativa.
  • executar e interpretar os testes estatísticos omnibus e post hoc.
  • executar e interpretar estatísticas de tamanho de efeito.

Preparação

Os exemplos aqui apresentados estão disponíveis. Caso queira usá-los, crie um projeto, coloque o arquivo desta aula e os seguintes arquivos na pasta do mesmo:

Sobre os métodos tradicionais

http://unusual-cars.com/wp-content/uploads/2016/01/Ford-Model-T-1908.jpg

ANOVA (sigla, do inglês, Analysis of Variance) utiliza a estatística \(F\).

  • A ANOVA de Fisher é uma extensão do teste \(t\), com:
    • VI (variável independente) nominal com três ou mais categorias (politômica nominal, fator)
    • VD (variável dependente) quantitativa (unifatorial, One way ANOVA)
  • Pode ser:
    • independente (fator entre participantes)
    • relacionada (fator intraparticipantes)

O tamanho do efeito é expresso pelo eta ao quadrado (\(\eta^2\)).


Tradicionalmente exigia-se, para a ANOVA clássica, que a distribuição populacional da variável fosse normal para cada condição e que existisse homocedasticidade entre os grupos. Estas premissas eram testadas:

  • Normalidade da VD (teste de Shapiro-Wilk) se a amostra fosse pequena.
  • Homocedasticidade (teste de Levene) se os grupos fossem desbalanceados.

Quando as premissas não eram atendidas, recomendava-se o teste não paramétrico H de Kruskal-Wallis.


ANOVA é relativamente robusta às violações das suposições de normalidade, homocedasticidade e outliers.

Se o tamanho da amostra é pequeno, ou há heterocedasticidade, desbalanceamento, ou assimetria da VD, deve-se considerar a execução de uma ANOVA de Welch, ou ANOVA não-paramétrica ou ANOVA com reamostragem (bootstrapping).

- é necessário seguir as premissas com todo o rigor?

Ilustramos com uma analogia, recordando o equilíbrio de Hardy-Weinberg.


The Hardy-Weinberg Law of Genetic Equilibrium

In 1908 G. Hardy and W. Weinberg independently proposed that the frequency of alleles and genotypes in a population will remain constant from generation to generation if the population is stable and in genetic equilibrium. Five conditions are required in order for a population to remain at Hardy-Weinberg equilibrium:

  • A large breeding population
  • Random mating
  • No change in allelic frequency due to mutation
  • No immigration or emigration
  • No natural selection


Violando as condições de H-W, em simulação com:

  • 500 indivíduos,
  • mutação alterando a frequência dos alelos,
  • selecionando para o alelo A.
source("HardyWeinberg.R")

Observamos, por exemplo: As condições matemáticas de população muito grande (ou infinita), sem mutação e sem seleção garantem a população no equilíbrio de Hardy-Weinberg; porém, dentro de certos limites, ainda encontramos predição aceitável para o Equilíbrio de Hardy-Weinberg e poderíamos utilizá-los em populações biológicas.

- o que é independência?

Um delineamento independente é aquele em que cada unidade experimental só é submetida a uma condição experimental.

Por exemplo, três grupos com 5 indivíduos cada um (total de 15 participantes, nomeados de A a J) são alocados assim:

Condição 1 Condição 2 Condição 3
A F K
B G L
C H M
D I N
E J O

Em um delineamento relacionado, o mesmo indivíduo é submetido a todas as condições experimentais. Os mesmos 15 participantes seriam alocados assim:

Condição 1 Condição 2 Condição 3
A A A
B B B
C C C
D D D
E E E
F F F
G G G
H H H
I I I
J J J
K K K
L L L
M M M
N N N
O O O

- o que é balanceamento?

Em delineamentos em que os grupos são balanceados, o número de participantes em cada condição é aproximadamente igual. Por exemplo, com 15 participantes em ANOVA independente, desbalanceada, poderíamos ter algo como.

Condição 1 Condição 2 Condição 3
A H K
B I L
C J M
D N
E O
F
G

Em um delineamento intra-participantes o desbalanceamento ocorre quando um dos participantes deixa de aparecer em uma das condições experimentais, por exemplo:

Condição 1 Condição 2 Condição 3
A A A
B B
C C C
D D D
E E E
F F F
G G G
H H H
I I I
J J
K K K
L L L
M M M
N N N
O O O

Basta uma participação faltando (missing) para caracterizar o desbalanceamento. O procedimento com as estatísticas padrão do R não funcionarão e temos que usar um modelo com efeitos aleatórios; lidaremos com esta situação adiante.

Distribuição \(F\) de Fisher-Snedecor

Esta distribuição é uma generalização da estatística \(t\) para três ou mais grupos. Descreve probabilidades que dependem da razão entre duas variâncias, e portanto considera graus de liberdade (\(\nu\), letra grega ni ) para o numerador e para o denominador.

Para a comparação de três ou mais condições, com base na estatística \(F\), os graus de liberdade do numerador dependem do número de condições e os do denominador dependem do tamanho da amostra.

Familiarize-se com a distribuição \(F\), observando Animacao_F.R:

  • \(F\) é um valor maior que zero.
  • sob \(H_0\) a distribuição \(F\) tem um parâmetro de não centralidade (ncp) igual a zero; sob \(H_1\) o parâmetro de não centralidade é maior do que zero.
  • esta distribuição é assimétrica.
  • consideramos somente a cauda superior.
  • localize \(\alpha\) e \(\beta\).
  • há dois valores para graus de liberdade: para o numerador (número de grupos - 1) e denominador (número de sujeitos - número de grupos).
  • observe o que acontece com a distribuição sob \(H_0\) e com o valor de \(F\) crítico (\(F_c\)) à medida que os graus de liberdade aumentam.

O delineamento dos estudos, o tipo de variável e, consequentemente, a estatística adequada mudam, mas o problema é sempre o mesmo: incerteza porque lidamos com uma amostra.

A diferença, aqui, é que teremos que lidar com 3 ou mais amostras simultaneamente.

A hipótese nula é pela igualdade de todos os \(m\) grupos. caso não rejeitemos \(H_0\):

Quando rejeitamos \(H_0\) assumimos que todas as condições são diferentes entre si:
ou

que pelo menos uma destoe das demais (i.e., basta que uma das condições ser diferente das demais para rejeitarmos \(H_0\).), por exemplo:

teste \(t\) (quando eram duas turmas)

Brendon Small e Coach McGuirk fazem com que seus alunos do SNAP-Ed mantenham diários do que comem por uma semana e depois calculem a ingestão diária de sódio em miligramas.

Desde que as classes receberam diferentes programas de educação nutricional, eles querem ver se a ingestão média de sódio é a mesma para as duas turmas.

Aplicamos o teste \(t\) de Satterthwaite (a.k.a. Welch), robusto e default do R:

# dados
library(readxl)
Dtfrm <- read_excel("Nutricao.xlsx")
# significancia estatistica
t_out <- t.test(Sodium ~ Instructor, data = Dtfrm)
print(t_out)

    Welch Two Sample t-test

data:  Sodium by Instructor
t = 0.76722, df = 34.893, p-value = 0.4481
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -67.91132 150.41132
sample estimates:
mean in group Brendon Small mean in group Coach McGuirk 
                    1287.50                     1246.25 

testes \(t\) quando temos três ou mais turmas

Suponha que uma terceira classe junte-se ao experimento (Melissa Robins). Podemos, então, verificar as diferenças destas três condições experimentais com testes \(t\), comparando-se Brendon Small com Coach McGuirk, Brendon Small com Melissa Robins e Coach McGuirk com Melissa Robins?

Quantos testes \(t\) precisaremos de acordo com o número de grupos?

É a combinatória do número de grupos (\(m\)) dois a dois: \[{m \choose 2} = { {m!} \over {2! (m-2)!} }\]

  • para 2 -grupos: \({2 \choose 2} = { {2!} \over {2! (2-2)!} } = 1~\text{teste}\)
  • para 3 grupos: \({3 \choose 2} = { {3!} \over {2! (3-2)!} } = 3~\text{testes}\)
  • para 4 grupos: \({4 \choose 2} = { {4!} \over {2! (4-2)!} } = 6~\text{testes}\)
  • para 6 grupos: \({6 \choose 2} = { {6!} \over {2! (6-2)!} } = 15~\text{testes}\)

Graficamente vemos:

# Testes_t_2a2.R
# Ilustra o numero de testes t necessarios
# para comparar grupos par a par

combinacoes <- c()
grupos <- c()
for (g in 2:20)
{
  grupos <- c(grupos,g)
  combinacoes <- c(combinacoes,choose(g,2))
}
plot(grupos,combinacoes,
     xlab="Número de grupos",
     ylab="Número de testes t feitos par a par",
     lwd=2,
     type="b")

Embora o número cresça rapidamente, podemos não nos impressionar pois temos computadores para o trabalho repetitivo. Há, porém, um problema mais grave: as probabilidades dos erros do tipo I (\(\alpha\)) e do tipo II (\(\beta\)) se acumulam. Quando escolhemos \(\alpha\), a probabilidade de rejeitar incorretamente a hipótese nula, o máximo valor que \(p\) pode assumir, a probabilidade deste erro para \(m\) grupos é:

\[P(\alpha | m) = 1-(1-\alpha)^{m \choose 2}\]

Também gostaríamos de manter a probabilidade de rejeitar corretamente a hipótese nula, o poder do teste, \(1-\beta\); para \(m\) grupos é:

\[P(\beta | m) = (1-\beta) ^ {m \choose 2}\]

Graficamente podemos observar o que acontece com o número crescente de pares de teste \(t\) necessários, considerando os tradicionais \(\alpha=0.05\) e \(\beta=0.1\):

# Testes_t_AlfaPoder.R
# Ilustra o valor de alfa e 1-beta
# com numero crescente de testes 

source("friendlycolor.R")

cor_alfa <- friendlycolor(8)
cor_poder <- friendlycolor(20)

alfa_escolhido <- 0.05
beta_escolhido <- 0.20
alfa <- c()
poder <- c()
grupos <- c()
for (g in 2:20)
{
  grupos <- c(grupos,g)
  combinacoes <- choose(g,2)
  
  alfa <- c(alfa,1-(1-alfa_escolhido)^combinacoes)
  poder <- c(poder,(1-beta_escolhido)^combinacoes)
}
plot(grupos, poder,
     xlab="Número de grupos",
     ylab="Probabilidade",
     ylim=c(0,1),
     lwd=2, lty=2, col=cor_poder,
     type="b")
lines(grupos,alfa,
      lwd=2, lty=1, col=cor_alfa, 
      type="b")
legend("right",
       c("alfa", "poder"),
       col=c(cor_alfa,cor_poder),
       lwd=c(2,2),
       lty=c(1,2),
       box.lwd=0, bg="transparent")

Portanto, a probabilidade de erro do tipo I cresce rapidamente para quase \(100\%\) e o poder de seus testes combinados vai para zero. Na prática, dificilmente teremos mais do que 6 grupos mas, ainda assim, teríamos \(\alpha \approx 53.7\%\) e poder \(1-\beta \approx 3.5\%\), valores totalmente inaceitáveis para uma boa análise.

Podemos resolver com vários testes \(t\)? A resposta é:

NÃO!

Precisamos testar tudo simultaneamente para
manter \(\alpha\) no valor pretendido e
preservar \(1-\beta\).

Análise da variância

Esta análise, que tem o acrônimo ANOVA (do inglês, Analysis of Variance) utiliza apenas variâncias, mas…

… as conclusões são sobre as médias dos \(m\) grupos.

Como?

Princípio do teste ANOVA unifatorial

Consideremos que existam três grupos com distribuição normal. De cada um deles, retiramos uma amostra: Não havendo problemas, as três amostras reproduzem a distribuição, média e variância das populações das quais se originaram.

A distribuição \(F\) é dada por \[F = {{S_E^2}\over{S_D^2}}\] onde \(S_E^2\) é a variância entre os grupos e \(S_D^2\) é a variância dentro dos grupos.

Variância entre (between) os grupos

A variância entre os grupos presume que as médias amostrais (\(\bar{X}_m\)) refletem as respectivas médias populacionais (\(\mu_m\)): Como sempre, não temos esta certeza e lidamos somente com a informação das amostras, \(\bar{X}_m\):

Cada média é um número e, portanto, podemos calcular a variância destes números: Sendo assim, para \(F = {{S_E^2}\over{S_D^2}}\), a estatística \(F\) aumenta quando a variância entre os grupos (ou condições) aumentar.

Variância dentro (within) dos grupos

A variância dentro dos grupos é uma medida da variância total, desconsiderando a média de cada uma das condições. Cada amostra tem sua própria distribuição (presumivelmente, reflexo da distribuição da população de onde veio): Como sempre, não temos esta certeza e lidaremos somente com a informação das amostras, \(\S_m\):

Mesclando as três distribuições, estimamos a variância dentro dos grupos, \(S_D^2\), uma medida de quanto, como um todo, a variável é dispersa:

Caso a variância em cada condição seja maior,

esta variância será refletida em \(S_D^2\): Sendo assim, para \(F = {{S_E^2}\over{S_D^2}}\), a estatística \(F\) diminui quando a variância dentro os grupos (ou condições) aumentar.

Comportamento de \(F = {{S_E^2}\over{S_D^2}}\)

É fácil imaginar o comportamento da estatística \(F\) combinando-se o que pode acontecer com \(S_E^2\) e \(S_D^2\):




Desta forma, Ronald Fisher inventou uma forma de comparar médias entre várias condições, simultaneamente, utilizando somente a comparação entre duas variâncias, com o numerador refletindo a dispersão das médias e o denominador refletindo a dispersão do fenômeno em estudo.

Métodos Robustos

https://performancedrive.com.au/icon-land-rover-defender-90-6-2-chev-v8-1606/

Os métodos aqui apresentados são robustos à falta de homocedasticidade. A ausência de normalidade é mais difícil de ser contornada em amostras pequenas. Podemos prescindir da suposição de normalidade para amostras maiores (tamanhos mínimos são vistos adiante) utilizando o teorema central do limite.

ANOVA para condições independentes

Utilizada quando os participantes são avaliados em somente uma das condições experimentais, i.e., um delineamento independente ou entre participantes.

- suposições

  • independência entre as unidades experimentais,
  • normalidade da VD em todas as condições,
  • homocedasticidade da VD entre todas as condições.

Normalidade e balanceamento:

“O teste t é ainda um teste válido, mesmo com modestas violações na suposição de normalidade, particularmente quando os tamanhos dos grupos são iguais e existe um número razoável de participantes em cada grupo; por “razoável” entendemos que, em um delineamento completamente entre participantes, deve haver pelo menos 12 participantes por grupo[…]” (Dancey & Reidy, 2019, p. 472)

Então: - se o número de unidades experimentais é igual a pelo menos 12 em cada grupo, a suposição de normalidade não é necessária.

  • se os grupos são balanceados, homocedasticidade não é necessária.

O balanceamento não precisa ser estrito:

“Consideram-se grupos de dimensão semelhante quando o quociente entre a maior dimensão e a menor for inferior a 1,5.” (Pestana, M. & Gageiro, J. (2008) Análise de dados para Ciências Sociais: a complementaridade do SPSS. Lisboa: Sílabo, p. 278)

Então:

  • Se a maior quantidade de unidades observacionais numa condição NÃO superar 1,5 vezes a condição de menor quantidade de unidades observacionais, então a suposição de homocedasticidade não precisa ser considerada para o teste ANOVA unifatorial independente.


Homocedasticidade:

Em geral, quando você tiver tamanhos de amostras iguais, a suposição de heterocedasticidade não será um grande problema.
- Dancey & Reidy (2019), p. 472-3
  • ANOVA unifatorial independente de Fisher para condições balanceadas é adequada na situação de heterocedasticidade populacional.
  • ANOVA unifatorial independente de Welch para condições desbalanceadas é adequada na situação de heterocedasticidade populacional.
A heterocedasticidade ocorre heuristicamente se (maior desvio-padrão)/(menor desvio-padrão) é maior que 2.
- Johnson, R. & Wichern, D. (2007) Applied Multivariate Statistical Analysis. 6 th ed. NJ: Prentice-Hall, p. 291
- Moore, D. (1995) The basic practice of statistics. New York: W. H. Freeman and Company,
citado por Norusis, M. (2009) PASW Statistics 18: Statistics Procedures Companion. NJ: Prentice-Hall, p. 148

- situação

O SNAP-Ed (Supplemental Nutrition Assistance Program Education) é um programa baseado em evidências que ajuda as pessoas a terem uma vida mais saudável.

O SNAP-Ed ensina às pessoas que usam ou qualificam para o SNAP uma boa nutrição e como fazer com que o seu dinheiro de alimentação se estenda ainda mais.

Os participantes do SNAP-Ed também aprendem a ser fisicamente ativos.

- planejamento

Brendon Small, Coach McGuirk e Melissa Robins fazem com que seus alunos do SNAP-Ed mantenham diários do que comem por uma semana e depois calculem a ingestão diária de sódio em miligramas.

Desde que as classes receberam diferentes programas de educação nutricional, eles querem ver se a ingestão média de sódio é a mesma para as três turmas.

- hipóteses nula e alternativa

As três classes receberam diferentes programas de educação nutricional. A ingestão média de sódio é a mesma (populacionalmente) para os três programas?

\[H_0: \mu_{\text{Small}} = \mu_{\text{McGuirk}} = \mu_{\text{Robins}}\] \[H_1: \exists |\mu_i \ne \mu_j; i \ne j; i,j=1,2,3\]


Esta é uma forma matemática para escrever a hipótese alternativa, lida como “Existe pelo menos alguma média \(\mu_i\), diferente de outra média \(\mu_j\), com \(i\) e \(j\) assumindo os valores \(1\), \(2\) ou \(3\)”. É dizer que \(\mu_1 \ne \mu_2\) OU \(\mu_1 \ne \mu_3\) OU \(\mu_2 \ne \mu_3\).

Aparece, por aí, como: \[H_1: \mu_1 \ne \mu_2 \ne \mu_3\] tentando indicar as diferenças, mas esta forma não funciona porque sugere que a rejeição de \(H_0\) implica em que todos os grupos sejam diferentes entre si, como \(\mu_1 \ne \mu_2\) E \(\mu_1 \ne \mu_3\) E \(\mu_2 \ne \mu_3\).

Basta que uma condição tenha média estatisticamente diferente das demais para que se rejeite \(H_0\) com ANOVA. Em português:

\[H_1: \text{Existe pelo menos um grupo diferente de algum outro.}\]

Operacionalização do teste

ANOVA unifatorial independente de Fisher com ajuste para heterocedasticidade de White

A estatística descritiva e inferencial pode ser obtida com ANOVA1f_indep_FisherWhite_sodio.R:

# ANOVA1f_indep_FisherWhite_sodio.R
# para ajustar este RScript para outros dados
# troque a planilha xlsx e substitua as palavras
# Instructor pela nova VI (fator)
# Sodium pela nova VD (resposta)

library(psych)
library(lattice)
library(car)
library(multcomp)
library(gplots)
library(emmeans)
library(readxl)

source("eiras_plotIC.R")

# suppress warnings
options(warn=-1)

TH <- read_excel("Nutricao3.xlsx")
TH$Instructor <- factor(TH$Instructor, levels=unique(TH$Instructor))

print(with(TH, psych::describeBy(Sodium,Instructor,digits=2)))
boxplot(Sodium~Instructor,data=TH,
        ylab=names(TH)[which(names(TH)=="Sodium")],
        xlab=names(TH)[which(names(TH)=="Instructor")]
)
print(grf <- lattice::xyplot(Sodium~Instructor, data=TH, type=c("p","a"),
                             jitter.x=TRUE, col="black"))
with(TH, gplots::plotmeans(Sodium~Instructor,
                           error.bars="conf.int", level=.95,
                           connect=FALSE,
                           ylab=names(TH)[which(names(TH)=="Sodium")],
                           xlab=names(TH)[which(names(TH)=="Instructor")],
                           main="IC95%",
                           barcol="black"))
car::densityPlot(Sodium~Instructor, data=TH, rug=TRUE, from=0, normalize=TRUE,
                 na.rm=TRUE, ylab="Densidade", col=c("black", "black", "black"))
cat("\n")

cat("\nANOVA unifatorial independente de Fisher",
    "\ncom ajuste para heterocedasticidade de White:\n\n")
alfa <- 0.05
VD <- names(TH)[which(names(TH)=="Sodium")]
VI <- names(TH)[which(names(TH)=="Instructor")]
cat("VD =", VD,"\n")
cat("Fator =", VI,"\n")
cat("\nAnalise de significancia estatistica: teste omnibus\n")
modelo <- lm(Sodium~Instructor, data=TH)
print(res <- car::Anova(modelo, type=2, white.adjust=TRUE))
cat("\nAnalise de significancia pratica: tamanho de efeito\n")
F <- res$F[1]
dfn <- res$Df[1]
dfd <- res$Df[2]
eta2 <- dfn*F/(dfn*F+dfd)
if (0 <= eta2 & eta2 < 0.1) {geta2 <- "minimo"}
if (0.1 <= eta2 & eta2 < 0.6) {geta2 <- "pequeno"}
if (0.6 <= eta2 & eta2 < 0.14) {geta2 <- "intermediario"}
if (0.14 <= eta2 & eta2 <= 1.0) {geta2 <- "grande"}
cat("- eta^2 =", eta2, "\nGrau", geta2,
    "de explicacao da variancia da VD", VD,"pela VI", VI,"\n")
f2 <- eta2/(1-eta2) # tamanho de efeito f de Cohen
ncp <- dfd*f2 # parametro de nao-centralidade
fc <- qf(1-alfa, dfn, dfd, 0)
p <- 1-pf(F,dfn,dfd,0)
if (p < 1e-4)
{
  p <- sprintf("%.2e",p)
} else
{
  p <- sprintf("%.4f",p)
}
f <- seq(0,2*ncp,0.01)
densf <- df(f, dfn, dfd, 0)
plot(f, densf, xlab="F", ylab="densidade", lwd=2, type="l")
densf <- df(f, dfn, dfd, ncp)
lines(f,densf, lwd=2, lty=2)
abline(v=fc, lty=3)
abline(v=F, lty=4)
legend("topright",
       c("H0", "Obs", 
         paste("Fc(",dfn,",",dfd,") = ",round(fc,3),sep=""), 
         paste("Fobs = ",round(F,3),"\n",
               "p = ",p,sep="") 
         ), 
       lwd=c(2,2,1,1), lty=c(1,2,3,4))
cat("\n\nSelecao de modelo\n")
R2aj <- (F-1)/((F-1)+dfd+1)
cat("- R^2 ajustado =", R2aj, "\n")
omega2 <- (F-1)/((F-1)+dfd+2)
cat("- omega^2 = ", omega2,"\n\n")
cat(paste("Teste post hoc\n"))
print(EMM <- emmeans::emmeans(modelo, "Instructor"))
print(grf <- plot(EMM, colors = "black",
                  main="Estimated Marginal Means",
                  xlab=VD,
                  ylab=VI))

# Testes post hoc
cat("\n\nTestes post hoc:\n")

# nomes do fator encurtados
TH$Instructor <- as.character(TH$Instructor)
fatores <- unique(as.character(TH$Instructor))
letra <- "A"
legenda <- c()
cat ("\nLegenda:\n")
for( f in 1:length(fatores))
{
  cat("\t",letra," ... ",fatores[f],"\n",sep="")
  legenda <- c(legenda,paste(letra," ... ",fatores[f],"\n",sep=""))
  TH$Instructor[TH$Instructor==fatores[f]] <- letra
  ascii <- strtoi(charToRaw(letra),16L)
  letra <- rawToChar(as.raw(ascii+1))
}
TH$Instructor <- as.factor(TH$Instructor)
modelo <- lm(Sodium~Instructor, data=TH)
mc.tukey <- multcomp::glht(modelo, linfct = mcp(Instructor = "Tukey"))
print(mcs.tukey <- summary(mc.tukey, test=adjusted("bonferroni")))
multcomp::cld(mcs.tukey, level=alfa, decreasing=TRUE)
# plot(mc.tukey,las=3)
t.ic <- mcs.tukey$test$qfunction(0.95)
df_plot <- data.frame(
  names (mcs.tukey$test$coefficient),
  as.vector(mcs.tukey$test$coefficients),
  as.vector(mcs.tukey$test$coefficients)-t.ic*as.vector(mcs.tukey$test$sigma),
  as.vector(mcs.tukey$test$coefficients)+t.ic*as.vector(mcs.tukey$test$sigma),
  as.vector(mcs.tukey$test$pvalues)
)
eiras_plotIC(df_plot,
             main="95% family-wise confidence level",
             xlab="Difference",
             usecolor="n"
             )
legend("topleft",legenda,lwd=0,lty=0,cex=0.6,box.lwd=0, border="transparent", bg="transparent")
mc.dunnett <- multcomp::glht(modelo, linfct = mcp(Instructor = "Dunnett"))
mcs.dunnett <- summary(mc.dunnett, test=adjusted("bonferroni"))
print(mcs.dunnett)
# plot(mc.dunnett,las=3)
t.ic <- mcs.dunnett$test$qfunction(0.95)
df_plot <- data.frame(
  names (mcs.dunnett$test$coefficient),
  as.vector(mcs.dunnett$test$coefficients),
  as.vector(mcs.dunnett$test$coefficients)-t.ic*as.vector(mcs.dunnett$test$sigma),
  as.vector(mcs.dunnett$test$coefficients)+t.ic*as.vector(mcs.dunnett$test$sigma),
  as.vector(mcs.dunnett$test$pvalues)
)
eiras_plotIC(df_plot,
             main="95% family-wise confidence level",
             xlab="Difference",
             usecolor="n"
)
legend("topleft",legenda,lwd=0,lty=0,cex=0.6,box.lwd=0, border="transparent", bg="transparent")

# enable warnings
options(warn=0)

Obtendo:

Loading required package: carData

Attaching package: 'car'
The following object is masked from 'package:psych':

    logit
Loading required package: mvtnorm
Loading required package: survival
Loading required package: TH.data
Loading required package: MASS

Attaching package: 'TH.data'
The following object is masked from 'package:MASS':

    geyser

Attaching package: 'gplots'
The following object is masked from 'package:stats':

    lowess

 Descriptive statistics by group 
group: Brendon Small
   vars  n   mean     sd median trimmed    mad min  max range skew
X1    1 20 1287.5 193.73   1300 1284.38 166.79 950 1700   750 0.12
   kurtosis    se
X1    -0.46 43.32
-------------------------------------------------------- 
group: Coach McGuirk
   vars  n    mean     sd median trimmed    mad  min  max range skew
X1    1 20 1246.25 142.41 1212.5 1240.62 148.26 1000 1525   525  0.3
   kurtosis    se
X1    -0.85 31.84
-------------------------------------------------------- 
group: Melissa Robins
   vars  n    mean     sd median trimmed    mad min  max range skew
X1    1 20 1123.75 143.15 1112.5 1120.31 185.32 900 1400   500 0.08
   kurtosis    se
X1     -1.1 32.01



ANOVA unifatorial independente de Fisher 
com ajuste para heterocedasticidade de White:

VD = Sodium 
Fator = Instructor 

Analise de significancia estatistica: teste omnibus
Coefficient covariances computed by hccm()

Analysis of Deviance Table (Type II tests)

Response: Sodium
           Df      F   Pr(>F)   
Instructor  2 5.5748 0.006148 **
Residuals  57                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Analise de significancia pratica: tamanho de efeito
- eta^2 = 0.1636056 
Grau grande de explicacao da variancia da VD Sodium pela VI Instructor 


Selecao de modelo
- R^2 ajustado = 0.0731098 
- omega^2 =  0.07195982 

Teste post hoc
 Instructor     emmean   SE df lower.CL upper.CL
 Brendon Small    1288 36.1 57     1215     1360
 Coach McGuirk    1246 36.1 57     1174     1319
 Melissa Robins   1124 36.1 57     1051     1196

Confidence level used: 0.95 
Registered S3 methods overwritten by 'ggplot2':
  method         from 
  [.quosures     rlang
  c.quosures     rlang
  print.quosures rlang



Testes post hoc:

Legenda:
    A ... Brendon Small
    B ... Coach McGuirk
    C ... Melissa Robins

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = Sodium ~ Instructor, data = TH)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
B - A == 0   -41.25      51.09  -0.807  1.00000   
C - A == 0  -163.75      51.09  -3.205  0.00664 **
C - B == 0  -122.50      51.09  -2.398  0.05939 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- bonferroni method)


     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: lm(formula = Sodium ~ Instructor, data = TH)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
B - A == 0   -41.25      51.09  -0.807  0.84559   
C - A == 0  -163.75      51.09  -3.205  0.00443 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- bonferroni method)

Analise a saída deste RScript e encontre:

  • a estatística descritiva (textual e gráficos)
  • a ANOVA unifatorial
    • encontre os graus de liberdade \(F\)(2, 57).
    • o valor de \(F\) calculado, igual a 5.5748, NA.
    • o valor de \(p\) correspondente, igual a 0.006148, NA.
    • o tamanho de efeito, \(\eta^2\)=0.1636056.
  • os testes post hoc com a legenda:
         A ... Brendon Small
 
         B ... Coach McGuirk
 
         C ... Melissa Robins
 
  • e seus resultados em:
    • contrastes de Tukey (textual e gráfico) com os respectivos valores \(p\): 1, 0.0066393, 0.0593866
    • contrastes de Dunnett (textual e gráfico) com os respectivos valores \(p\): 0.8455883, 0.0044262
    • note que \(p=1\) não pode existir, é um arredondamento; por este motivo, no gráfico, fizemos uma correção.

Existem muitas famílias de testes para análise post hoc disponíveis.

Uma família de testes para análise post-hoc mais conservadora é adequada quando há um grande número de comparações como, por exemplo, o das diferenças honestamente significantes (DHS) de Tukey.

“O Teste proposto por Tukey (1953) é também conhecido como teste de Tukey da diferença honestamente significativa (honestly significant difference)(HSD) e teste de Tukey da diferença totalmente significativa (wholly significant difference)(WSD). É um teste exato em que, para a família de todas as \(c=\frac{1}{2}k(k-1)\) comparações duas a duas, a taxa de erro da família dos testes (FWER) é exatamente \(\alpha\) (e o intervalo de confiança é exatamente \(1-\alpha\)). Métodos de comparações múltiplas exatos são raros. O teste de Tukey tem sido mostrado analiticamente ótimo, no sentido que, entre todos os procedimentos que resultam em intervalos de confiança com mesmo tamanho para todas diferenças duas a duas com coeficiente de confiança da família de pelo menos \(1-\alpha\), o teste de Tukey resulta em intervalos menores. Isso quer dizer que, se a família consiste em todas comparações duas a duas e o teste de Tukey pode ser usado, ele resultará em intervalos menores que qualquer outro método de comparação múltipla de uma etapa.”

Exemplo de relatório para a ANOVA unifatorial independente

Há 60 participantes no estudo balanceado com três condições independentes, sendo que nenhuma delas é de controle. Os participantes foram distribuídos aleatoriamente e balanceadamente nos três grupos. As suposições de normalidade e homocedasticidade da VD foram consideradas válidas. As médias amostrais brutas mostram que menor ingesta de sódio foi observada entre estudantes submetidos ao programa aplicado por Melissa Robins. A ingesta de sódio dos estudantes de Brendon Small e do Coach McGuirk são semelhantes. A análise de variância de um fator fixo entre participantes mostrou que o efeito fixo Sódio é estatisticamente significante, pois o teste omnibus produziu F(2;57) = 5,57 e p = 0,00615. O tamanho do efeito do fator Instrutor é expresso por eta ao quadrado, sendo que seu valor é igual a 0,16. Portanto, 16% da variância da ingesta de sódio é explicada pelo programa adotado pelos respectivos instrutores. Os testes post hoc de Tukey-HSD confirmaram que as diferenças entre os programas adotados por Robins e Small, e entre os adotados por Robins e McGuirk são estatisticamente significantes. Não se observou diferença estatisticamente significante entre os programas adotados por Small e McGuirk.

ANOVA unifatorial independente de Welch

O RScript ANOVA1f_indep_Welch_sodio.R é uma modificação mais robusta do anterior.

As linhas de código que produzem a ANOVA de Fisher-White:

modelo <- lm(Sodium~Instructor, data=TH)
print(res <- car::Anova(modelo, type=2, white.adjust=TRUE))

são substituídas por:

print(res <- jmv::anovaOneW(formula = Sodium~Instructor, data=TH,
                       desc=TRUE, descPlot = FALSE, phMethod ='gamesHowell',
                       phMeanDif = TRUE, phTest=TRUE, phFlag=TRUE))

O objeto res, que recebe o que retorna da função jmv::anovaOneW() é diferente de car::Anova() e, portanto, na sequência do RScript, a captação dos valores foram devidamente ajustados. Esta função da library jmv já executa os testes post hoc.

Executando:

# ANOVA1f_indep_Welch_sodio.R
# para ajustar este RScript para outros dados
# troque a planilha xlsx e substitua as palavras
# Instructor pela nova VI (fator)
# Sodium pela nova VD (resposta)

library(psych)
library(lattice)
library(car)
library(gplots)
library(jmv)
library(readxl)

# suppress warnings
options(warn=-1)

TH <- read_excel("Nutricao3.xlsx")
TH$Instructor <- factor(TH$Instructor, levels=unique(TH$Instructor))

print(with(TH, psych::describeBy(Sodium,Instructor,digits=2)))
boxplot(Sodium~Instructor,data=TH,
        ylab=names(TH)[which(names(TH)=="Sodium")],
        xlab=names(TH)[which(names(TH)=="Instructor")]
)
print(grf <- lattice::xyplot(Sodium~Instructor, data=TH, type=c("p","a"),
                             jitter.x=TRUE, col="black"))
with(TH, gplots::plotmeans(Sodium~Instructor,
                           error.bars="conf.int", level=.95,
                           connect=FALSE,
                           ylab=names(TH)[which(names(TH)=="Sodium")],
                           xlab=names(TH)[which(names(TH)=="Instructor")],
                           main="IC95%",
                           barcol="black"))
car::densityPlot(Sodium~Instructor, data=TH, rug=TRUE, from=0, normalize=TRUE,
                 na.rm=TRUE, ylab="Densidade", col=c("black", "black", "black"))
cat("\n")

cat("\nANOVA unifatorial independente de Welch\n\n")
alfa <- 0.05
VD <- names(TH)[which(names(TH)=="Sodium")]
VI <- names(TH)[which(names(TH)=="Instructor")]
cat("VD =", VD,"\n")
cat("Fator =", VI,"\n")
cat("\nAnalise de significancia estatistica: teste omnibus\n")
print(res <- jmv::anovaOneW(formula = Sodium~Instructor, data=TH,
                       desc=TRUE, descPlot = FALSE, phMethod ='gamesHowell',
                       phMeanDif = TRUE, phTest=TRUE, phFlag=TRUE))
cat("\nAnalise de significancia pratica: tamanho de efeito\n")
F <- as.numeric(res$anova$asDF[2])
dfn <- as.numeric(res$anova$asDF[3])
dfd <- as.numeric(res$anova$asDF[4])
eta2 <- dfn*F/(dfn*F+dfd)
if (0 <= eta2 & eta2 < 0.1) {geta2 <- "minimo"}
if (0.1 <= eta2 & eta2 < 0.6) {geta2 <- "pequeno"}
if (0.6 <= eta2 & eta2 < 0.14) {geta2 <- "intermediario"}
if (0.14 <= eta2 & eta2 <= 1.0) {geta2 <- "grande"}
cat("- eta^2 =", eta2, "\nGrau", geta2,
    "de explicacao da variancia da VD", VD,"pela VI", VI,"\n")
f2 <- eta2/(1-eta2) # tamanho de efeito f de Cohen
ncp <- dfd*f2 # parametro de nao-centralidade
fc <- qf(1-alfa, dfn, dfd, 0)
p <- 1-pf(F,dfn,dfd,0)
if (p < 1e-4)
{
  p <- sprintf("%.2e",p)
} else
{
  p <- sprintf("%.4f",p)
}
f <- seq(0,2*ncp,0.01)
densf <- df(f, dfn, dfd, 0)
plot(f, densf, xlab="F", ylab="densidade", lwd=2, type="l")
densf <- df(f, dfn, dfd, ncp)
lines(f,densf, lwd=2, lty=2)
abline(v=fc, lty=3)
abline(v=F, lty=4)
legend("topright",
       c("H0", "Obs", 
         paste("Fc(",dfn,",",round(dfd,3),") = ",round(fc,3),sep=""), 
         paste("Fobs = ",round(F,3),"\n",
               "p = ",p,sep="") 
         ), 
       lwd=c(2,2,1,1), lty=c(1,2,3,4))
cat("\n\nSelecao de modelo\n")
R2aj <- (F-1)/((F-1)+dfd+1)
cat("- R^2 ajustado =", R2aj, "\n")
omega2 <- (F-1)/((F-1)+dfd+2)
cat("- omega^2 = ", omega2,"\n\n")
# enable warnings
options(warn=0)

Obtém-se:


Attaching package: 'jmv'
The following object is masked from 'package:psych':

    pca
The following object is masked from 'package:stats':

    anova

 Descriptive statistics by group 
group: Brendon Small
   vars  n   mean     sd median trimmed    mad min  max range skew
X1    1 20 1287.5 193.73   1300 1284.38 166.79 950 1700   750 0.12
   kurtosis    se
X1    -0.46 43.32
-------------------------------------------------------- 
group: Coach McGuirk
   vars  n    mean     sd median trimmed    mad  min  max range skew
X1    1 20 1246.25 142.41 1212.5 1240.62 148.26 1000 1525   525  0.3
   kurtosis    se
X1    -0.85 31.84
-------------------------------------------------------- 
group: Melissa Robins
   vars  n    mean     sd median trimmed    mad min  max range skew
X1    1 20 1123.75 143.15 1112.5 1120.31 185.32 900 1400   500 0.08
   kurtosis    se
X1     -1.1 32.01



ANOVA unifatorial independente de Welch

VD = Sodium 
Fator = Instructor 

Analise de significancia estatistica:

     One-Way ANOVA (Welch's)
            dep F[welch] df1[welch] df2[welch]    p[welch]
"Sodium" Sodium 5.765465          2   37.39637 0.006569226

    Group Descriptives
                          dep          group num    mean       sd       se
"SodiumBrendon Small"  Sodium  Brendon Small  20 1287.50 193.7341 43.32026
"SodiumCoach McGuirk"  Sodium  Coach McGuirk  20 1246.25 142.4123 31.84435
"SodiumMelissa Robins" Sodium Melissa Robins  20 1123.75 143.1495 32.00920

    Games-Howell Post-Hoc Test – Sodium
                                                                          
                                Brendon Small Coach McGuirk Melissa Robins
                                                                          
 Brendon Small  Mean difference -             41.25         163.75        
                t-value         -             0.7672        3.0401        
                df              -             34.8931       34.9827       
                p-value         -             0.7253        0.0121        
                                                                          
 Coach McGuirk  Mean difference -             -             122.5         
                t-value         -             -             2.7131        
                df              -             -             37.999        
                p-value         -             -             0.0263        
                                                                          
 Melissa Robins Mean difference -             -             -             
                t-value         -             -             -             
                df              -             -             -             
                p-value         -             -             -             
                                                                          

Analise de significancia pratica: tamanho de efeito
- eta^2 = 0.2356748 
Grau grande de explicacao da variancia da VD Sodium pela VI Instructor 



Selecao de modelo
- R^2 ajustado = 0.1104092 
- omega^2 =  0.1079091 

Analise a saída deste RScript e encontre:

  • a mesma estatística descritiva do RScript anterior.
  • os valores da ANOVA de Welch, com o valor \(F\text{ calculado}\)=(2,37.4)=5.7655 e o valor-\(p\)=0.0065692 (note o número fracionário para os graus de liberdade do denominados, da mesma forma que acontece com o teste \(t\) de Welch).
  • os testes post hoc, na forma de uma matriz relacionando os pares de programas (Small, McGuirk e Robins), mostrando:
    • a diferença média
    • o valor da estatística \(t\) (são feitos testes \(t\) par-a-par)
    • os graus de liberdade fracionários (\(t\) de Welch)
    • os valores \(p\), já corrigidos para serem interpretados em comparação com \(\alpha\)=0.05.
  • o valor estimado para o tamanho de efeito do teste omnibus, dado por \(\eta^2\)=0.2356748.
  • note que este procedimento não fornece os IC95% pós-modelo; pretendemos incorporar isto depois

A conclusão, neste exemplo, é similar à obtida com o teste de ANOVA de Fisher-White: o programa adotado por Melissa Robins obteve ingestas de sódio significantemente menores que a dos outros dois instrutores.

ANOVA unifatorial por bootstrapping

O RScript ANOVA1f_indep_Bootstrap_sodio.R é a versão com bootstrapping do pacote lmboot.

O RScript que implementa ANOVA desta maneira é:

# ANOVA1f_indep_Bootstrap_sodio.R
# para ajustar este RScript para outros dados
# troque a planilha xlsx e substitua as palavras
# Instructor pela nova VI (fator)
# Sodium pela nova VD (resposta)

library(psych)
library(lattice)
library(car)
library(gplots)
library(lmboot)
library(readxl)

# suppress warnings
options(warn=-1)

TH <- read_excel("Nutricao3.xlsx")
TH$Instructor <- factor(TH$Instructor, levels=unique(TH$Instructor))

print(with(TH, psych::describeBy(Sodium,Instructor,digits=2)))
boxplot(Sodium~Instructor,data=TH,
        ylab=names(TH)[which(names(TH)=="Sodium")],
        xlab=names(TH)[which(names(TH)=="Instructor")]
)
print(grf <- lattice::xyplot(Sodium~Instructor, data=TH, type=c("p","a"),
                             jitter.x=TRUE, col="black"))
with(TH, gplots::plotmeans(Sodium~Instructor,
                           error.bars="conf.int", level=.95,
                           connect=FALSE,
                           ylab=names(TH)[which(names(TH)=="Sodium")],
                           xlab=names(TH)[which(names(TH)=="Instructor")],
                           main="IC95%",
                           barcol="black"))
car::densityPlot(Sodium~Instructor, data=TH, rug=TRUE, from=0, normalize=TRUE,
                 na.rm=TRUE, ylab="Densidade", col=c("black", "black", "black"))
cat("\n")

cat("\nANOVA unifatorial por bootstrapping\n\n")
alfa <- 0.05
VD <- names(TH)[which(names(TH)=="Sodium")]
VI <- names(TH)[which(names(TH)=="Instructor")]
cat("VD =", VD,"\n")
cat("Fator =", VI,"\n")
cat("\nAnalise de significancia estatistica: teste omnibus\n")

# remova seed para ter resultados variados
bootsamples <- 1e6
modelo_boot <- lmboot::ANOVA.boot(Sodium~Instructor, 
                                  B = bootsamples, type = "residual", 
                                  wild.dist = "normal",  
                                  # seed = 123,
                                  data = TH, 
                                  keep.boot.resp = FALSE)
cat(paste("F(",dfn,",",round(dfd,3),") = ",round(F,5),", p = ",round(modelo_boot$`p-values`,5),"\n", sep=""))
cat(paste("(",bootsamples," bootstrap samples)\n", sep=""))
cat("\nAnalise de significancia pratica: tamanho de efeito\n")
F <- qf(1-modelo_boot$`p-values`, modelo_boot$df[1], modelo_boot$df[2])
dfn <- as.numeric(modelo_boot$df[1])
dfd <- as.numeric(modelo_boot$df[2])
eta2 <- dfn*F/(dfn*F+dfd)
if (0 <= eta2 & eta2 < 0.1) {geta2 <- "minimo"}
if (0.1 <= eta2 & eta2 < 0.6) {geta2 <- "pequeno"}
if (0.6 <= eta2 & eta2 < 0.14) {geta2 <- "intermediario"}
if (0.14 <= eta2 & eta2 <= 1.0) {geta2 <- "grande"}
cat("- eta^2 =", eta2, "\nGrau", geta2,
    "de explicacao da variancia da VD", VD,"pela VI", VI,"\n")
f2 <- eta2/(1-eta2) # tamanho de efeito f de Cohen
ncp <- dfd*f2 # parametro de nao-centralidade
fc <- qf(1-alfa, dfn, dfd, 0)
p <- modelo_boot$`p-values`
if (p < 1e-4)
{
  p <- sprintf("%.2e",p)
} else
{
  p <- sprintf("%.4f",p)
}
f <- seq(0,2*ncp,0.01)
densf <- df(f, dfn, dfd, 0)
plot(f, densf, xlab="F", ylab="densidade", lwd=2, type="l")
densf <- df(f, dfn, dfd, ncp)
lines(f,densf, lwd=2, lty=2)
abline(v=fc, lty=3)
abline(v=F, lty=4)
legend("topright",
       c("H0", "Obs", 
         paste("Fc(",dfn,",",round(dfd,3),") = ",round(fc,3),sep=""), 
         paste("Fobs = ",round(F,3),"\n",
               "p = ",p,sep="") 
         ), 
       lwd=c(2,2,1,1), lty=c(1,2,3,4))
cat("\n\nSelecao de modelo\n")
R2aj <- (F-1)/((F-1)+dfd+1)
cat("- R^2 ajustado =", R2aj, "\n")
omega2 <- (F-1)/((F-1)+dfd+2)
cat("- omega^2 = ", omega2,"\n\n")
# enable warnings
options(warn=0)

Obtendo-se:


 Descriptive statistics by group 
group: Brendon Small
   vars  n   mean     sd median trimmed    mad min  max range skew
X1    1 20 1287.5 193.73   1300 1284.38 166.79 950 1700   750 0.12
   kurtosis    se
X1    -0.46 43.32
-------------------------------------------------------- 
group: Coach McGuirk
   vars  n    mean     sd median trimmed    mad  min  max range skew
X1    1 20 1246.25 142.41 1212.5 1240.62 148.26 1000 1525   525  0.3
   kurtosis    se
X1    -0.85 31.84
-------------------------------------------------------- 
group: Melissa Robins
   vars  n    mean     sd median trimmed    mad min  max range skew
X1    1 20 1123.75 143.15 1112.5 1120.31 185.32 900 1400   500 0.08
   kurtosis    se
X1     -1.1 32.01



ANOVA unifatorial por bootstrapping

VD = Sodium 
Fator = Instructor 

Analise de significancia estatistica: teste omnibus
F(2,37.396) = 5.76546, p = 0.00635
(1e+06 bootstrap samples)

Analise de significancia pratica: tamanho de efeito
- eta^2 = 0.1626598 
Grau grande de explicacao da variancia da VD Sodium pela VI Instructor 



Selecao de modelo
- R^2 ajustado = 0.07253934 
- omega^2 =  0.07139764 

Verifique na saída:

  • compare com as duas soluções anteriores.
  • note que este procedimento não fornece os IC95% pós-modelo; pretendemos incorporar isto depois

ANOVA independente, sem os dados brutos

É muito comum, em publicações, que somente tenhamos acesso às medidas-resumo (número de participantes, média, desvio-padrão e correlação). Nestes casos, os RScripts acima não são utilizáveis.

Para fazer os testes quando os dados brutos não estão disponíveis, criamos os seguintes scripts:

ANOVA unifatorial relacionada ou para medidas repetidas

Utilizada quando os participantes são avaliados sob todas as condições experimentais, i.e., um delineamento intraparticipantes.

Neste delineamento a variação entre os grupos não é devida às diferenças individuais (cada participante em cada grupo é o mesmo), e a fórmula para o cálculo da estatística \(F\) leva este fato em consideração.

O delineamento com medidas repetidas leva a um teste com maior poder.

No entanto, como o mesmo participante é controle de si mesmo e submetido a todas as condições experimentais, sugere-se que esta ordem de submissão, se possível, seja aleatorizada.

- suposições

– as diferenças dos valores das VDs são independentes entre as unidades observacionais. – as diferenças dos valores das VDs têm distribuição normal multivariada. – esfericidade: homocedasticidade das variâncias das diferenças das VDs (uma explicação sem muita matemática está em https://en.wikipedia.org/wiki/Mauchly%27s_sphericity_test)

“O teste t é ainda um teste válido, mesmo com modestas violações na suposição de normalidade, particularmente quando os tamanhos dos grupos são iguais e existe um número razoável de participantes em cada grupo; por “razoável” entendemos que [para um delineamento] completamente intraparticipantes, [deve haver] pelo menos 22 participantes ao todo.” (Dancey & Reidy, 2019, p. 472)

- situação


Usaremos os mesmos dados do exemplo anterior, mas imaginando que as 20 medidas de ingesta de sódio sejam do mesmo participante, submetido aos três diferentes programas educacionais. Então, a normalidade multivariada das três diferenças não pode ser assumida, pois há 20 unidades observacionais no estudo.

O SNAP-Ed (Supplemental Nutrition Assistance Program Education) é um programa baseado em evidências que ajuda as pessoas a terem uma vida mais saudável.

O SNAP-Ed ensina às pessoas que usam ou qualificam para o SNAP uma boa nutrição e como fazer com que o seu dinheiro de alimentação se estenda ainda mais.

Os participantes do SNAP-Ed também aprendem a ser fisicamente ativos.

- planejamento

Brendon Small, Coach McGuirk e Melissa Robins fazem com que seus alunos do SNAP-Ed mantenham diários do que comem por uma semana e depois calculem a ingestão diária de sódio em miligramas.

Estudantes atenderam os diferentes programas de educação nutricional suscessivamente, e os instrutores querem ver se a ingestão média de sódio é a mesma quando cada um dos três programas foi seguido.

- hipóteses nula e alternativa

O delineamento do estudo é diferente, mas as hipóteses são as mesmas:

\[H_0: \mu_{\text{Small}} = \mu_{\text{McGuirk}} = \mu_{\text{Robins}}\] \[H_1: \exists |\mu_i \ne \mu_j; i \ne j; i,j=1,2,3\]

Operacionalização do teste

ANOVA unifatorial relacionada ou para medidas repetidas, balanceadas

A planilha a ser utilizada é Nutricao3par.xlsx. Verifique a coluna Student, que foi alterada para indicar o participante nos três programas (compare com Nutricao3.xlsx, usada para ANOVA independente).

A libraries necessárias são lmerTest e ez e o trecho que aplica o método estatístico é:

res <- lmerTest::lmer(Sodium~Instructor + (1|Student), data=TH, REML=TRUE)

Executando:

# ANOVA1f_dep_balanc_sodio.R
# para ajustar este RScript para outros dados
# troque a planilha xlsx e substitua as palavras
# Instructor pela nova VI (fator)
# Sodium pela nova VD (resposta)

library(psych)
library(lattice)
library(car)
library(lmerTest)
library(ez)
library(readxl)

# suppress warnings
options(warn=-1)

TH <- read_excel("Nutricao3par.xlsx")
TH$Instructor <- factor(TH$Instructor, levels=unique(TH$Instructor))

print(with(TH, psych::describeBy(Sodium,Instructor,digits=2)))
boxplot(Sodium~Instructor,data=TH,
        ylab=names(TH)[which(names(TH)=="Sodium")],
        xlab=names(TH)[which(names(TH)=="Instructor")]
)
print(grf <- lattice::xyplot(Sodium~Instructor, data=TH, type=c("p","a"),
                             jitter.x=TRUE, col="black"))
car::densityPlot(Sodium~Instructor, data=TH, rug=TRUE, from=0, normalize=TRUE,
                 na.rm=TRUE, ylab="Densidade", col=c("black", "black", "black"))
cat("\n")

cat(paste("\nTeste omnibus por GLMM com efeito aleatorio\n"))
alfa <- 0.05
VD <- names(TH)[which(names(TH)=="Sodium")]
VI <- names(TH)[which(names(TH)=="Instructor")]
cat("VD =", VD,"\n")
cat("Fator =", VI,"\n")
cat("\nAnalise de significancia estatistica: teste omnibus\n")
res <- lmerTest::lmer(Sodium~Instructor + (1|Student), data=TH, REML=TRUE)
cat("\n"); print(res1 <- stats::anova(res)); cat("\n")
print(lmerTest::rand(res))
cat(paste("\nTeste post hoc\n\n"))
print(out <- lmerTest::difflsmeans(res))
dfn <- as.numeric(res1[3])
dfd <- as.numeric(res1[4])
F <- as.numeric(res1[5])
cat("\nAnalise de significancia pratica: tamanho de efeito\n")
eta2 <- dfn*F/(dfn*F+dfd)
if (0 <= eta2 & eta2 < 0.1) {geta2 <- "minimo"}
if (0.1 <= eta2 & eta2 < 0.6) {geta2 <- "pequeno"}
if (0.6 <= eta2 & eta2 < 0.14) {geta2 <- "intermediario"}
if (0.14 <= eta2 & eta2 <= 1.0) {geta2 <- "grande"}
cat("\neta^2 =", eta2, "\nGrau", geta2,
    "de explicacao da variancia da VD", VD,"pela VI", VI,"\n")
f2 <- eta2/(1-eta2) # tamanho de efeito f de Cohen
ncp <- dfd*f2 # parametro de nao-centralidade
fc <- qf(1-alfa, dfn, dfd, 0)
p <- 1-pf(F,dfn,dfd,0)
if (p < 1e-4)
{
  p <- sprintf("%.2e",p)
} else
{
  p <- sprintf("%.4f",p)
}
f <- seq(0,2*ncp,0.01)
densf <- df(f, dfn, dfd, 0)
plot(f, densf, xlab="F", ylab="densidade", lwd=2, type="l")
densf <- df(f, dfn, dfd, ncp)
lines(f,densf, lwd=2, lty=2)
abline(v=fc, lty=3)
abline(v=F, lty=4)
legend("topright",
       c("H0", "Obs", 
         paste("Fc(",dfn,",",round(dfd,3),") = ",round(fc,3),sep=""), 
         paste("Fobs = ",round(F,3),"\n",
               "p = ",p,sep="") 
         ), 
       lwd=c(2,2,1,1), lty=c(1,2,3,4))
cat("\n\nSelecao de modelo\n")
R2aj <- (F-1)/((F-1)+dfd+1)
cat("- R^2 ajustado =", R2aj, "\n")
omega2 <- (F-1)/((F-1)+dfd+2)
cat("- omega^2 = ", omega2,"\n\n")
cat("Teste omnibus por GLM univariado\n")
# ezANOVA: apenas para medidas repetidas balanceadas
print(res <- ez::ezANOVA(data=TH, dv=Sodium, wid=Student, within=Instructor,
                         detailed=TRUE, type=3))
# enable warnings
options(warn=0)

Obtém-se:

Loading required package: lme4
Loading required package: Matrix
Registered S3 methods overwritten by 'lme4':
  method                          from
  cooks.distance.influence.merMod car 
  influence.merMod                car 
  dfbeta.influence.merMod         car 
  dfbetas.influence.merMod        car 

Attaching package: 'lmerTest'
The following object is masked from 'package:lme4':

    lmer
The following object is masked from 'package:stats':

    step

 Descriptive statistics by group 
group: Brendon Small
   vars  n   mean     sd median trimmed    mad min  max range skew
X1    1 20 1287.5 193.73   1300 1284.38 166.79 950 1700   750 0.12
   kurtosis    se
X1    -0.46 43.32
-------------------------------------------------------- 
group: Coach McGuirk
   vars  n    mean     sd median trimmed    mad  min  max range skew
X1    1 20 1246.25 142.41 1212.5 1240.62 148.26 1000 1525   525  0.3
   kurtosis    se
X1    -0.85 31.84
-------------------------------------------------------- 
group: Melissa Robins
   vars  n    mean     sd median trimmed    mad min  max range skew
X1    1 20 1123.75 143.15 1112.5 1120.31 185.32 900 1400   500 0.08
   kurtosis    se
X1     -1.1 32.01



Teste omnibus por GLMM com efeito aleatorio
VD = Sodium 
Fator = Instructor 

Analise de significancia estatistica: teste omnibus

Type III Analysis of Variance Table with Satterthwaite's method
           Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
Instructor 290146  145073     2    38  30.581 1.217e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA-like table for random-effects: Single term deletions

Model:
Sodium ~ Instructor + (1 | Student)
              npar  logLik    AIC    LRT Df Pr(>Chisq)    
<none>           5 -352.02 714.05                         
(1 | Student)    4 -375.21 758.42 46.376  1  9.763e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Teste post hoc

Least Squares Means table:

                                                   Estimate Std. Error df
InstructorBrendon Small - InstructorCoach McGuirk   41.2500    21.7807 38
InstructorBrendon Small - InstructorMelissa Robins 163.7500    21.7807 38
InstructorCoach McGuirk - InstructorMelissa Robins 122.5000    21.7807 38
                                                   t value    lower
InstructorBrendon Small - InstructorCoach McGuirk   1.8939  -2.8426
InstructorBrendon Small - InstructorMelissa Robins  7.5181 119.6574
InstructorCoach McGuirk - InstructorMelissa Robins  5.6243  78.4074
                                                      upper  Pr(>|t|)    
InstructorBrendon Small - InstructorCoach McGuirk   85.3426   0.06587 .  
InstructorBrendon Small - InstructorMelissa Robins 207.8426 4.949e-09 ***
InstructorCoach McGuirk - InstructorMelissa Robins 166.5926 1.866e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Confidence level: 95%
  Degrees of freedom method: Satterthwaite 

Analise de significancia pratica: tamanho de efeito

eta^2 = 0.6167848 
Grau grande de explicacao da variancia da VD Sodium pela VI Instructor 



Selecao de modelo
- R^2 ajustado = 0.4313252 
- omega^2 =  0.4251262 

Teste omnibus por GLM univariado
$ANOVA
       Effect DFn DFd        SSn       SSd          F            p p<.05
1 (Intercept)   1  19 89182041.7 1307541.7 1295.91189 6.002751e-19     *
2  Instructor   2  38   290145.8  180270.8   30.58049 1.217342e-08     *
        ges
1 0.9835909
2 0.1631905

$`Mauchly's Test for Sphericity`
      Effect         W            p p<.05
2 Instructor 0.3889714 0.0002038255     *

$`Sphericity Corrections`
      Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF]
2 Instructor 0.6207214 3.881854e-06         * 0.6427565 2.772133e-06
  p[HF]<.05
2         *

Analise a saída deste RScript e encontre:

  • a estatística descritiva (sempre a mesma).
  • observe o valor de \(F\text{ calculado}\) e valor-\(p\) correspondentes (neste exemplo fictício, os dados são os mesmos da ANOVA independente).
  • localize e interprete os testes post hoc; compare com os valores-\(p\) da ANOVA independente.
  • aparece um novo trecho, “Teste omnibus por GLM univariado”.

ANOVA unifatorial relacionada ou para medidas repetidas, desbalanceadas

Vamos usar os mesmos dados, mas eliminar uma única medida (removemos o estudante a de Melissa Robins, planilha Nutricao3parD.xlsx).

Para lidar com esta situação, o RScript execute:

# ANOVA1f_dep_desbalanc_sodio.R
# para ajustar este RScript para outros dados
# troque a planilha xlsx e substitua as palavras
# Instructor pela nova VI (fator)
# Sodium pela nova VD (resposta)

library(psych)
library(lattice)
library(car)
library(lmerTest)
library(readxl)

# suppress warnings
options(warn=-1)

TH <- read_excel("Nutricao3parD.xlsx")
TH$Instructor <- factor(TH$Instructor, levels=unique(TH$Instructor))

print(with(TH, psych::describeBy(Sodium,Instructor,digits=2)))
boxplot(Sodium~Instructor,data=TH,
        ylab=names(TH)[which(names(TH)=="Sodium")],
        xlab=names(TH)[which(names(TH)=="Instructor")]
)
print(grf <- lattice::xyplot(Sodium~Instructor, data=TH, type=c("p","a"),
                             jitter.x=TRUE, col="black"))
car::densityPlot(Sodium~Instructor, data=TH, rug=TRUE, from=0, normalize=TRUE,
                 na.rm=TRUE, ylab="Densidade", col=c("black", "black", "black"))
cat("\n")

cat(paste("\nTeste omnibus por GLMM com efeito aleatorio\n"))
alfa <- 0.05
VD <- names(TH)[which(names(TH)=="Sodium")]
VI <- names(TH)[which(names(TH)=="Instructor")]
cat("VD =", VD,"\n")
cat("Fator =", VI,"\n")
cat("\nAnalise de significancia estatistica: teste omnibus\n")
res <- lmerTest::lmer(Sodium~Instructor + (1|Student), data=TH, REML=TRUE)
cat("\n"); print(res1 <- stats::anova(res)); cat("\n")
print(lmerTest::rand(res))
cat(paste("\nTeste post hoc\n\n"))
print(out <- lmerTest::difflsmeans(res))
dfn <- as.numeric(res1[3])
dfd <- as.numeric(res1[4])
F <- as.numeric(res1[5])
cat("\nAnalise de significancia pratica: tamanho de efeito\n")
eta2 <- dfn*F/(dfn*F+dfd)
if (0 <= eta2 & eta2 < 0.1) {geta2 <- "minimo"}
if (0.1 <= eta2 & eta2 < 0.6) {geta2 <- "pequeno"}
if (0.6 <= eta2 & eta2 < 0.14) {geta2 <- "intermediario"}
if (0.14 <= eta2 & eta2 <= 1.0) {geta2 <- "grande"}
cat("\neta^2 =", eta2, "\nGrau", geta2,
    "de explicacao da variancia da VD", VD,"pela VI", VI,"\n")
f2 <- eta2/(1-eta2) # tamanho de efeito f de Cohen
ncp <- dfd*f2 # parametro de nao-centralidade
fc <- qf(1-alfa, dfn, dfd, 0)
p <- 1-pf(F,dfn,dfd,0)
if (p < 1e-4)
{
  p <- sprintf("%.2e",p)
} else
{
  p <- sprintf("%.4f",p)
}
f <- seq(0,2*ncp,0.01)
densf <- df(f, dfn, dfd, 0)
plot(f, densf, xlab="F", ylab="densidade", lwd=2, type="l")
densf <- df(f, dfn, dfd, ncp)
lines(f,densf, lwd=2, lty=2)
abline(v=fc, lty=3)
abline(v=F, lty=4)
legend("topright",
       c("H0", "Obs", 
         paste("Fc(",dfn,",",round(dfd,3),") = ",round(fc,3),sep=""), 
         paste("Fobs = ",round(F,3),"\n",
               "p = ",p,sep="") 
         ), 
       lwd=c(2,2,1,1), lty=c(1,2,3,4))
cat("\n\nSelecao de modelo\n")
R2aj <- (F-1)/((F-1)+dfd+1)
cat("- R^2 ajustado =", R2aj, "\n")
omega2 <- (F-1)/((F-1)+dfd+2)
cat("- omega^2 = ", omega2,"\n\n")
# enable warnings
options(warn=0)

É o mesmo RScript usado para a ANOVA balanceada, exceto pelas linhas finais (“Teste omnibus por GLM univariado”). Obtém-se:


 Descriptive statistics by group 
group: Brendon Small
   vars  n   mean     sd median trimmed    mad min  max range skew
X1    1 20 1287.5 193.73   1300 1284.38 166.79 950 1700   750 0.12
   kurtosis    se
X1    -0.46 43.32
-------------------------------------------------------- 
group: Coach McGuirk
   vars  n    mean     sd median trimmed    mad  min  max range skew
X1    1 20 1246.25 142.41 1212.5 1240.62 148.26 1000 1525   525  0.3
   kurtosis    se
X1    -0.85 31.84
-------------------------------------------------------- 
group: Melissa Robins
   vars  n    mean     sd median trimmed    mad min  max range skew
X1    1 19 1135.53 136.76   1125 1132.35 148.26 925 1400   475 0.08
   kurtosis    se
X1    -1.04 31.37



Teste omnibus por GLMM com efeito aleatorio
VD = Sodium 
Fator = Instructor 

Analise de significancia estatistica: teste omnibus

Type III Analysis of Variance Table with Satterthwaite's method
           Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
Instructor 254746  127373     2 37.086  27.369 4.998e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA-like table for random-effects: Single term deletions

Model:
Sodium ~ Instructor + (1 | Student)
              npar  logLik    AIC    LRT Df Pr(>Chisq)    
<none>           5 -345.58 701.16                         
(1 | Student)    4 -368.17 744.34 45.177  1    1.8e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Teste post hoc

Least Squares Means table:

                                                   Estimate Std. Error
InstructorBrendon Small - InstructorCoach McGuirk   41.2500    21.5730
InstructorBrendon Small - InstructorMelissa Robins 157.5073    21.9807
InstructorCoach McGuirk - InstructorMelissa Robins 116.2573    21.9807
                                                     df t value    lower
InstructorBrendon Small - InstructorCoach McGuirk  37.0  1.9121  -2.4601
InstructorBrendon Small - InstructorMelissa Robins 37.1  7.1657 112.9749
InstructorCoach McGuirk - InstructorMelissa Robins 37.1  5.2891  71.7249
                                                      upper  Pr(>|t|)    
InstructorBrendon Small - InstructorCoach McGuirk   84.9601   0.06362 .  
InstructorBrendon Small - InstructorMelissa Robins 202.0397 1.677e-08 ***
InstructorCoach McGuirk - InstructorMelissa Robins 160.7897 5.709e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Confidence level: 95%
  Degrees of freedom method: Satterthwaite 

Analise de significancia pratica: tamanho de efeito

eta^2 = 0.5961188 
Grau grande de explicacao da variancia da VD Sodium pela VI Instructor 



Selecao de modelo
- R^2 ajustado = 0.4091071 
- omega^2 =  0.4028569 

Exemplo de relatório para a ANOVA unifatorial relacionada

Vinte participantes foram selecionados para um estudo com delineamento intrapartipantes com três condições experimentais. A ordem de aplicação das três condições experimentais em cada participante foi aleatorizada. As médias amostrais brutas mostram que menor ingesta de sódio foi observada entre estudantes submetidos ao programa aplicado por Melissa Robins. A ingesta de sódio dos estudantes de Brendon Small e do Coach McGuirk são semelhantes. A ANOVA unifatorial para medidas repetidas mostrou que qualquer diferença entre as condições experimentais é improvável de ter ocorrido apenas por erro amostral, considerando a hipótese nula verdadeira, pois a estatística de teste multivariado observada é F(2;38) = 30,581 e o valor-p associado é igual a \(1.22 \cdot 10^{-8}\). O tamanho do efeito do fator intraparticipantes álcool é estimado pelo eta ao quadrado cujo valor indica que 61,7% da variância da ingesta de sódio é explicada pelo efeito do fator fixo Instrutor. A análise post-hoc confirmou que as diferenças entre entre os programas adotados por Robins e Small, e entre os adotados por Robins e McGuirk são estatisticamente significantes. Não se observou diferença estatisticamente significante entre os programas adotados por Small e McGuirk.